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r^ • ABSTRACT 

O 

O ■ To model the interaction between the solar wind and the interstellar wind, 

^ ■ magnetic fields must be included. Recently Opher et al. 2003 found that, by 

including the solar magnetic field in a 3D high resolution simulation using the 
University of Michigan BATS-R-US code, a jet-sheet structure forms beyond the 



c/3 . solar wind Termination Shock. Here we present an even higher resolution three- 

c^ ■ 

dimensional case where the jet extends for 150AU beyond the Termination Shock. 



We discuss the formation of the jet due to a de Laval nozzle effect and it's sub- 
sequent large period oscillation due to magnetohydrodynamic instabilities. To 
verify the source of the instability, we also perform a simplified two dimensional- 
geometry magnetohydrodynamic calculation of a plane fiuid jet embedded in a 
neutral sheet with the profiles taken from our 3D simulation. We find remarkable 
agreement with the full three-dimensional evolution. We compare both simula- 
tions and the temporal evolution of the jet showing that the sinuous mode is the 
dominant mode that develops into a velocity-shear-instability with a growth rate 
of 5 X 10~^sec~^ = 0.027 years'^. As a result, the outer edge of the heliosphere 
presents remarkable dynamics, such as turbulent fiows caused by the motion of 
the jet. Further study, e.g., including neutrals and the tilt of the solar rotation 
from the magnetic axis, is required before we can definitively address how this 
outer boundary behaves. Already, however, we can say that the magnetic field 
effects are a major player in this region changing our previous notion of how the 
solar system ends. 
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Subject headings: instabilities - interplanetary medium - ISM:kinematics and 
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1. Introduction 

One of the most exciting phenomena in astrophysical media are jets. Usually, they are 
collimated structures extending for several parsecs near black holes and active galaxies (e.g., 
Sahayanathan et al. 2003). Recently, we showed (Opher et al. 2003) that due the compression 
of the solar magnetic field, a jet-sheet structure forms in the region beyond the Termination 
Shock, where the solar wind becomes subsonic. This structure is narrow and collimated in 
the meridional plane (jet), but extends in the equatorial plane as a disk for approximately 
60° around the upstream direction. This result was obtained using the state-of-the art three 
dimensional magnetohydrodynamic code (BATS-R-US) with an unprecedented resolution to 
investigate the region where the solar wind meets the interstellar medium. It was motivated 
by the approach of the Voyager spacecraft to this region. Voyager is one of the most successful 
missions of NASA, traveling for more than 25 years with two spacecraft, respectively with a 
velocity of 3.5AU/year (35° above the ecliptic plane) and 3.1AU/year (48° below the ecliptic 
plane) . 

In this section we review the current state of modeling the interaction between the solar 
and the interstellar wind and discuss why magnetic effects are crucial, especially in the region 
between the Termination Shock and the Heliopause, the Heliosheath. 

Because the solar system travels in the interstellar medium with a velocity of approx- 
imately 25km /s (Frisch 1996), the solar wind collides with a supersonic interstellar wind. 
The basic structures that are formed by this interaction are: the Termination Shock (TS), 
the Heliopause (HP), and, possibly a Bow Shock (BS). The Termination Shock is the bound- 
ary where the supersonic solar wind becomes subsonic slowing down approaching the He- 
liopause. The Heliopause is the boundary separating the subsonic solar outflow downstream 
of the Termination Shock and the interstellar plasma flowing around it. The region inside 
the Heliopause is called the Heliosphere. If the flow of the interstellar wind is also super- 
alfvenic, there is a BS further out, where the interstellar flow becomes subsonic to avoid 
the heliospheric obstacle. In previous hydrodynamic models, the region beyond the Termi- 
nation Shock has constant plasma pressure and temperature and the heliospheric boundary 
was a smooth, rounded surface. This is not true if the solar magnetic field is included: the 
plasma pressure, temperature and density downstream the shock are not uniform and con- 
stant, and the heliospheric boundary is highly distorted from the rounded appearence of the 
hydrodynamic models. 
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The region between the Termination Shock and the Hehopause, the Hehosheath, is one of 
the most mysterious and unknown regions. The total plasma beta (the ratio between the sum 
of the thermal and ram pressure and the magnetic pressure) decreases as one passes through 
the shock; it is clear that the magnetic field will play a major role in that region (Before the 
TS, Ptotai ~ 65. After the TS Ptotai ~ 8. At the magnetic ridges jStotai < 1 ~ 0.1). The solar 
magnetic field reverses polarity at the heliospheric current sheet (HCS). If we neglect the 
tilt between the solar magnetic and rotation axis, the HCS remains in the equatorial plane. 
One of the major questions is how the HCS behaves beyond the Termination Shock, or in 
another words, how the current system in the Heliosphere closes. 

The interaction between the solar and interstellar wind is inherently a three dimensional 
problem because of the magnetic fields. Ideally, modeling this complex region properly 
would require inclusion of several ingredients: a) the ionized component from the solar and 
interstellar wind; b) the neutral atoms (coupled to the plasma by charge-exchange collisions); 
c) the magnetic field (of both the solar wind and the interstellar wind); d) the pickup ions; 
and e) cosmic rays (both galactic and anomalous components). 

For the interstellar wind, we know with certainty neither the magnitude nor the direction 
of the magnetic field. In the case of the neutral component, the associated particle large 
mean free path dictates a kinetic treatment. In addition, the diverse length scales present in 
the problem must be treated consistently (i.e. the numerical resolution of a simulation code 
must be capable of widely differing resolutions according to the different regions). The scales 
that we are interested in are on the order of tenths of an AU (for MHD instabilities and the 
heliospheric current sheet, for example) although the dimensions of the Heliosphere are on 
the order of a thousand AU. Therefore, another aspect that is important when modeling this 
interaction is the numerical resolution. 

There have been several numerical approaches to tackle this complicated interaction. 
Much work has focused on the careful treatment of the neutral component. Baranov & 
Malama (1993, 1995) (see also Baranov, Izmodenov & Malama (1998)) developed a kinetic 
description for the neutrals using a Monte Carlo approach. These studies, however, neglected 
magnetic fields. More recently, Baranov & Zaitsev (1995) and Izmodenov, Gruntman & 
Malama (2001) used the same description and included the interstellar magnetic field in a two 
dimensional model. Another alternative approach, still in two dimensions and neglecting the 
magnetic field, was developed by Pauls & Zank (1996) (see also (Zank et al. 1996; Pauls et al. 
1995)) who treated the neutrals as an additional fluid species. More recently, Muller, Zank & 
Lipatov (2000) used a particle-in-mesh approach as a kinetic approximation for the neutrals. 
Liewer et al. (1996) and Karmesin et al. (1995) used a gasdynamic approach treating the 
neutrals as a single fluid, also in a two dimensional study. The works by Ratkiewicz et al. 
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(1998) (see also (Ratkiewicz, Barnes, & Spreiter 2000; Ratkiewicz & McKenzie 2003)) and 
(Ratkiewicz & Ben-Jaffel 2002) were three-dimensional studies assuming that the neutrals 
have constant density, velocity and temperature. They also neglected the solar magnetic field 
and studied the effect of the interstellar magnetic field for several magnitudes and directions. 

The works by Pogorolev & Matsuda (2000), McNutt, Lyon & Goodrich (1998), Washimi 
& Tanaka (1996, 2001) and Linde et al. (1998) were the only ones that included both the 
solar and interstellar magnetic fields in a self consistent way in three-dimensional treatments. 
The neutrals were included as a single fiuid. The major drawback of these studies are: the 
inclusion of both of the interstellar and solar magnetic fields and the neutrals making it 
difficult to differentiate the effect of each component; and the limited spatial resolution. 

In our recent study (Opher et al. 2003) we performed a three dimensional study including 
the solar magnetic field. In order to study how the HCS behaves beyond the TS, we did 
not include the interstellar magnetic field to eliminate the effect of reconnection between the 
solar and interstellar fields and to isolate the effect of the solar magnetic field. We used an 
adaptive mesh refinement allowing us to get to spatial resolutions previously not obtained 
of the order of 1.5 AU at the HCS. We showed that a jet-sheet structure forms at the edge 
of the solar system. This jet-sheet oscillates up and down due to a velocity shear instability. 
In that study, we assume as a first approximation that the solar magnetic and the rotation 
axis are aligned. We also only treat the ionized component and neglect the effect of neutrals. 
Both effects, i.e., the inclusion of neutrals and the tilt of the dipole, need to be included to 
assess the overall structure and dynamics of the edge of our solar system. Spatial resolution 
was a key factor for resolving the jet-sheet structure at the edge of the solar system. We 
showed that resolutions of the order of 3 AU are insufficient; the jet at the current sheet is 
broadened and the current sheet remains in the equatorial plane as in the previous studies 
of Linde et al. (1998); Washimi & Tanaka (1996). 

Figure 1 shows the pressure contours of a meridional cut of our simulation with and 
without the inclusion of the solar magnetic field. The coordinate system is such that the 
interstellar wind is fiowing from the negative x-direction and the solar rotation axis is in 
the positive z-direction. The run presented in both cases is the coarse grid case (see Opher 
et al. (2003) and Appendix) with resolution of ~ 3 AU at the HCS. It is clear that in the 
Heliosheath the pressure is not uniform anymore, being enhanced in the region of the current 
sheet. Also, the HP is more pointed and the distance between the Heliopause and the Bow 
Shock increases (the HP location moves from 240 AU to 280 AU and the BS location moves 
from 358 AU to 440 AU), due to the presence of the jet. 

In the present study we revisit the work of Opher et al. (2003) in order to study the 
formation of the jet and the magnetohydrodynamic (MHD) instabilities driven by the jet- 
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sheet. We present a run where we extended the resolution to 0.75 AU at the HCS. (In the 
previous study we used a resolution of 1.5 AU at the current sheet, 3-4 times more refined 
than previous studies.) With such high resolution and extending the refined region, we were 
able to resolve the jet extending to 150 AU beyond the Termination Shock. We discuss on 
the effect of the de Laval nozzle on the formation of the jet and the width of the jet. 

The outline of the paper is the following: In the first section we describe briefly the 
model (for more details, see Opher et al. (2003)). In the second section, we present our 
major results for the run used and discuss the effect of the de Laval nozzle on formation of 
the jet. In the third section we discuss the width of the jet. The fourth section is dedicated 
to the discussion of the role of magnetohydrodynamic instabilities and comparison of the 
3D and 2D MHD code. The fifth section is dedicated to the conclusions and discussions. 
Finally, in the appendix we discuss the dependence of the width of the jet on resolution, 
showing that at the high resolution case, we obtained grid converged numerical solution. 



2. Description of the Model 

The code that we are using, BATS-R-US (Block Adaptive Tree Solar wind Roe type 
Upwind Scheme), was developed at the University of Michigan (Powell et al. 1999). The BATS- 
R-US code has been designed to capitalize on advances on massive parallel computers; the 
evolution of adaptive mesh refinement (AMR) techniques and advances in basic numerical 
methods particularly for hyperbolic conservations laws. The BATS-R-US MHD algorithm 
uses an upwind methods, approximate Riemann solvers, and limited solution reconstruction 
(see Powell et al. (1999)). It is a three dimensional magnetohydrodynamic code and the 
governing equations are the ideal magnetohydrodynamic equations. The resulting finite- 
volume scheme solves for the hydrodynamic and electromagnetic effects in a tightly coupled 
matter. It uses block-based solution for the adaptation on a cartesian mesh and is able to 
have more that 20 levels of refinement. It is a parallel multi-processor code with capability 
of running in different platforms, such as Beowulf (type PC clusters) and SGI origin 2000 
machines. 

In this study, as in Opher et al. (2003), we only treat the ionized components neglecting 
the neutrals. In order to eliminate the effect of reconnection at the heliopause, we did not 
include an interstellar magnetic field. The inner boundary was chosen at 30 AU and the 
outer boundary at 3450 x 4500 x 4500 AU. The physical boundary conditions are: For the 
solar wind (at 30 AU): Upiasma = 7.8 x lO'^ cm'^, T = 1.6 x 10^ K, v = 450 km/s. The 
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solar magnetic field at the inner boundary was taken as the Parker spiral magnetic field, 

— * 

B = sign{cosQ) 



-DO Cr — -Do 



r / \ r I u 



where i?o is taken as the inner boundary 30 AU , Ugw is the solar wind speed, 450 km/s, O is 
the polar angle of the field line, and Qq is the equatorial angular velocity of the Sun. At the 
equatorial plane, B = Bq = 2fiG. We chose the field polarity {sign{cosQ)) to match the 1996 
solar minimum. For the interstellar wind: Upiasma = O.OTcm'^, T = 10^ K, v = 25 km/s. 



2.1. Jet formation and properties: magnetic ridges and the de Laval nozzle 

In the present study we use higher spatial resolutions than the previously reported 
run (Opher et al. 2003) obtaining a resolution up to 0.75 AU at the current sheet. In the 
Appendix we comment on the effect of spatial resolution on the width of the jet. In this 
section, we briefiy summarize the major results of this run indicating the major features and 
differences from the previous case (with resolution of 1.5 AU at the HCS). Figure 2a shows 
the contours of the magnetic field at t = 1.3 x 10^ sec (scale ranging from — 0.3 riT). 
The coordinate system is the same as in Figure 1. Beyond the the Termination Shock (TS), 
the azimuthal magnetic field is compressed forming magnetic ridges (strong field, see region 
in red). The spatial resolution at the magnetic ridges is 3 AU, double of the resolution of 
previous works (Opher et al. 2003) and the magnetic field intensity ranges from — 0.4 uT 
(an increase from 0.24 uT from the previous study). Figure 2b shows the contours of the 
velocity Ux{iT^ the meridional plane) between the TS and the HP. In both, the black lines, 
are the fiow streamlines. 

The velocity of the jet is ~ 200 km/s extending for 150 AU beyond the TS, to x = 
—320 AU , in the direction of the BS (located at a: = —380 AU). Due to the higher resolution 
the turbulent vortices are more pronounced (with sizes of ~ 20 AU) and the structure is 
more complex than in Opher et al. (2003). The region where the turbulent vortices form has 
a resolution of 1.5 — 3.0 AU , while the one in Opher et al. (2003) has a resolution of 6.0 AU 
in the same region (see Appendix for details). The jet is unstable as in Opher et al. (2003). 
We discuss details of the instability in Section 4. 

The jet-sheet structure forms in the region of minimum of magnetic pressure. Figure 
2c shows that the fiow is in pressure balance after the shock. It portrays the vertical line 
cut downstream of the TS (at a; = 161 AU), at y = for thermal (dashed line), magnetic 
(dash- dot- dash line) and total (solid line) pressure. Downstream of the shock, where the 
flow decelerates, conservation of magnetic flux outside of the equatorial plane causes the 
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field to increase its magnitude further, while in the current sheet there is no such effect. As 
a resuh, the increased magnetic field above and below the ecliptic planes effectively pinches 
the sheet just beyond the Termination Shock, causing the stream lines in the subsonic regions 
to converge slightly. As pointed out by Washimi & Tanaka (2001), the enhanced magnetic 
ridges also create an obstacle in the post-shock flow, decelerating and deviating it to the 
flanks before reaching the Heliopause. This can be seen in Figure 2d where a line cut is 
taken, 60° above the ecliptic plane. The ram pressure (the green curve) reaches a minimum 
at the maximum intensity of the magnetic pressure, the magnetic ridges (red curve). (See 
also streamlines in Figure 2a). 

Figure 3 shows the ratio of velocity of the jet u with respect to the sound speed a at 
t = 1.3 X lO^sec, when the current sheet is still in the equatorial plane (green curve). For 
comparison, two other cases are presented: the refined (red) and coarse (blue) resolution 
cases. For the coarse case, the de Laval nozzle effect is negligble due to a broader jet. As we 
increase the resolution (red curve), the jet extends farther away, in the direction of the Bow 
Shock. In the refined case (as in Opher et al. (2003)) the jet loses it's power at a; = —230 AU 
(indicated by the green arrow) while the change in resolution (from 1.5 to 3.0 AU) occurs 
at = —320 AU . In the extended case the jet loses it's power at x = —320 AU (while the BS 
is at a; = —380 AU , indicated by the red arrow). The change in resolution (0.75 to 1.5 AU) 
occurs at —350 AU and, from 1.5 to 3.0 AU at —365 AU . This indicates that resolution 
plays a major role and the jet loses it's power due to numerical dissipation. 

The converging flow lines near the equatorial plane create a de Laval nozzle: a subsonic 
flow must accelerate where the streamlines converge and decelerate where flow lines diverge, 
while the opposite is true for a supersonic flow. Figure 4 shows a schematic drawing of a 
nozzle along with the streamlines where the jet form. The red arrow indicate the region 
where acceleration occurs. For a steady gas flow the Euler's equation can be written as 
(Landau & Lifshitz 1987) 

dv _ dp 

dr dr 

where F is the force. Consider a flow of gas through a nozzle with a variable cross section 

A[r) (as pictured in Figure 4a) Considering the pressure p = pRT, where R is the gas 

constant and the mass flux constant pvA{r) = constant, the equation above can be written 

as 

v^ - c,2 dv T IdA 



dr p A dr 



(3) 



For a narrowing nozzle, dA/dr < 0, and the flow velocity increases as long as it is 
subsonic f^ < 05^, reaching a maximum where the cross section is minimum. If the nozzle 



then widens again, the flow decelerates thereafter. However, if the flow accelerates sufficiently 
to become supersonic at the neck, then it will continue to accelerate where the nozzel widens 
{dA/dr > 0). A nozzle that ffist narrows and then widens again is called the de Laval 
nozzle. This is analogous to what happens in Parker's model of the solar wind at the critical 
point, where the effective nozzle, created by the combined effects of mass flux conservation 
and gravity, has a minimum cross section. In the jet formed here, the maximum speed just 
barely exceeds the sound speed, and decelerates thereafter, much as a breeze in solar wind 
theory. 

At the TS the flow velocity decreases to Mach number 0.55. Due to the acceleration of 
the flow past the TS, the Mach number increases to 1.1. The velocity of the jet is ~ 150 km/ s 
and remains almost constant for the extension of the jet (see Figure 3, where the velocity of 
the jet decreases only at a; = —320 AU due to the change in resolution). Fig. 4 shows the 
location of the de Laval nozzle effect. 



3. Width of the Jet 

What determines the scaling of the width of the jet? Does it depend on the grid 
resolution at the current sheet? Observations from Pioneer and Voyager indicate that the 
width of the HCS, at 1 AU is approximately 10, 000 km = 6.7x 10^^ AU (while a surrounding 
plasma sheet is thicker by a factor of 30 (Smith 2001; Winterhalter et al. 1994)) In terms 
of proton gyroradius the current sheet thickness is about 200 Rl (Winterhalter et al. 1994). 
For all cases, this is smaller than our grid resolution. 

How does the width of the current sheet scales with distance? Neglecting dynamic effects 
(such as turbulence) and kinetic effects such as drift instabilities, and assuming pressure 
balance, the width is determined by a balance between the thermal (inside the sheet) and 
the magnetic pressure outside. Combining the adiabatic energy and the energy equation, 

{u ■ V)p + 7p(V ■ m) = . (4) 

This equation can be written as 

dp ^ 

uq-t- + 7PV ■ M = , (5) 

where uq is the radial solar wind velocity, and r the radial distance from the Sun. Let the 
width of the current sheet be a, from Eq. (5), 

1 dp ^ 7 d{ar) 
p dr ar dr 



-9- 

that gives the dependence for the thermal pressure with distance as 

p = (ar)'"' , (7) 

where for an adiabatic gas 7 = 5/3. For a spherical expansion a oc r and the thermal pressure 
falls as r^^^^^. However, since the azimuthal magnetic field falls with distance as r~^, the 
magnetic pressure falls with r~^. Therefore, from the pressure equilibrium that must exist 
between the thermal pressure at the current sheet and the magnetic pressure outside the 
sheet, the dependence of the width with distance must be, 

a{r) oc r^/^ . (8) 

Using this scaling, for the obserevd width of 6.7 x 10^^ AU at 1 AU, the current sheet width 
at 150 AU is 2 X 10~^ AU which is much smaller than the grid resolution of our most resolved 
case. 

If the width of the HCS is determined by numerical diffusion, we can estimate the width 
as follows: For a steady state solution, the induction equation gives the width of the current 
sheet to be 

V -{uB- Bu) = rimV^B (9) 

or approximately 

Vm^^U,^^, (10) 

giving the approximately the width as 



^ , (11) 



U, 



i.e., the width of the HCS, a, in steady state, is determined by numerical diffusion rjm and 
the converging flows to the current sheet. The numerical diffusion 77^ will be proportional 
to the grid spatial resolution. If numerical diffusion determines width, the width decreases 
as resolution increases. However, this is not observed. 

Figure 5a shows the dependence of the jet width on the grid resolution in the region 
of the current sheet (—10 AU < 2; < 10 AU). Three cases are included: coarse, refined 
(same as Opher et al. 2003) and the super-refined-extended (the one presented at section 
2.1). For the coarse grid, the resolution at the current sheet is 6 AU , for the refined grid, 
1.5 AU and for the super-refined-extended grid, 0.75 AU . The blue and the red curves are, 
respectively, the half width of the velocity and the magnetic field, measured after the TS, 
at a; = —210 AU . The half widths of both the velocity and the magnetic field are measured 
half way of the maximum. For the three cases, the velocity (jet) profile follows closely the 
magnetic field (the current sheet) profile, being always slightly broader. As we increase the 



-10- 



resolution at the flanks of the jet, the two profiles approaches each other, almost overlapping. 
(For more details see Appendix). As we increase the resolution, the width of the velocity (jet) 
approaches the width of the magnetic field (current sheet). As we increase the resolution, 
the magnetic field (current sheet) width converges to a value of 4 AU. 

Figure 5b presents the the velocity profiles for the three cases: the coarse, the refined and 
the super-refined-extended cases. The red, green and blue curves, are respectively, the super- 
refined-extended, refined and coarse cases. The jet is much more resolved at the fianks as 
well as at the region of the current sheet in the super-refined-extended and the refined cases, 
compared to the coarse case. The profiles for the velocity for the super-refined-extended and 
the refined case are almost identical. We can see that the width of the jet, in the super- 
refined-extended case, is independent of the grid resolution and is determined by the physical 
conditions, rather than by numerical resolution or numerical diffusion. 



4. Instabilities - Comparison of 3D with 2D MHD code 

There have been several studies examining magnetized flows in which both sheared 
flows and sheared magnetic flelds are present (Einaudi 1999; Dahlburg, Keppens & Einaudi 
2001; Einaudi et al. 2001). A particular conflguration studied, motivated by the study of 
the heliospheric current sheet, was of a fluid jet embedded in a neutral sheet. Studies like 
Dahlburg, Boncinelli & Einaudi (1998); Einaudi et al. (2001) analyzed the behavior of a 
plane magnetized jet. Both study were two-dimensional cases (the justification came from 
laminar nature of the initial system and of the fact that the fastest growing linearly unstable 
modes are 2D). The configuration studied was a fluid jet in the z direction and the spatial 
coordinate along which the mean flow varies as y: Wo{y) = sech{y)ez embedded in a neutral 
sheet Bq = tanh{6y)ez, where Wq and Bq are the basic velocity fleld and magnetic fleld and 
6 is the magnetic shear width. These studies focus on the response of the wake-neutral or 
jet-neutral sheet to infinitesimal perturbations. The two unstable modes that exist in the 
fiuid limit (no magnetic field) are the varicose (sausage-like) and sinuous (kink-like) modes. 
These modes differ in the symmetry of the perturbed velocity component in the direction 
perpendicular to the flow. For the varicose mode, the perturbed velocity component is 
antisymmetric in y and, for the sinuous mode it is symmetric. In a magnetically dominated 
limit, only a tearing mode exists. This mode has the same symmetry as the varicose fluid 
mode. In the flow-dominated regime the sinuous mode has the fastest linear growth rate 
(faster than the varicose mode). Both modes are Kelvin- Helmholtz modes of a jet. The 
sinuous and varicose modes become magnetohydrodynamic modes for flnite values of the 
Alfven number A, the ratio of the characteristic Alfven speed to the characteristic flow 
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speed. 

In the 3D case presented in Section 2.1, v^ <^ lOO/cm/s and Vfiow ~ 200 km/s so 
the Alfven number, A <^ 1. At the location of the jet, there is no magnetic field and we 
expect the predominant mode to be the sinuous mode. The jet in the 3D case presented in 
Section 2.1 is unstable and oscillates up and down as in the case presented in Opher et al. 
(2003), although due to the higher resolution, much more structure can be seen (see attached 
material for a movie of Figure 2a). Figure 6a presents a measure of the displacement of the 
current sheet as a function of time. The red curve is the displacement in the middle of the 
jet, at a; = —225 AU, and the green curve the displacement at the nose, at a; = —290 AU. 
Indexes (I-III) indicate the different regions of time scales of oscillation. The displacement at 
the two locations of the jet are anti-correlated, implying that the wavelength of the oscillation 
A is < 70 AU. Examining Fig. 6b it is clear that the oscillation is a sinuous mode and that 
the wavelength is A ~ 23 AU. There are two distinct time scales. At the beginning of the 
displacement, a time scale of oscillation of ti = 2 x lO^sec (region I) and then t2 = 2x lO^sec 
(region II). Region III marks the exponential growth of the velocity-shear-Kelvin-Helmholtz 
instability. The growth rate, measured for curve III is F = 5.4x 10~^sec~^. (As a comparison, 
the growth rate measured for the refined case in Opher et al. (2003) was 1.6 x lO^^sec"^.) 

The Kelvin Helmholtz (KH) growth rate, for an infinitely thin layer is (Chandrasekhar 
1961) F = 0.5x \ k ■ Uq \ [1 — (2c^A; ■ Bq) /(k ■ Uq) Y^"^ ., where Uq is the velocity jump across 
the shear layer, ca is the Alfven speed, k is the direction of the wavevector (here x) and 
Bq is the direction of the magnetic field. Maximum growth occurs for k ~ 1/a, where a 
is the distance between the maximum velocity and zero. For the 3D run performed here, 
at a; = —190 AU , Uq ^ 150 km/s, and a ^ 20 AU. With these values the predicted KH 
growth rate, for planar geometry is 5 x 10^^ sec^^. 

To verify that this instability seen in the 3D MHD code is indeed caused by a KH type 
instability arising in the jet-sheet, we studied a much simpler 2D configuration, analogous to 
that studied by Einaudi (1999). In our case, the radial magnetic field in the fiow direction is 
negligible, while there is a strong sheared magnetic field perpendicular to this direction. The 
sheet thickness is much smaller than both the sheet length and it's longitudinal extension. 
Therefore, in the simplified geometry (Bettarini 2003), we model this by considering an initial 
equilibrium structure in cartesian geometry with velocity Vy that varies in the x direction 
and is uniform in the y and z directions. The jet-sheet is therefore infinite in the y — z 
plane. The line profiles for the density, temperature and magnetic field were extracted from 
the case described in Section 2.1. Figure 7a and 7b presents the initial profiles taken for 
the velocity of the jet VOy as a function of latitude x and for the magnetic field BOz as a 
function of the latitude x. 
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In particular, while in the 3D jet the temperature varies by a factor of almost 10, the 
density within the jet decreases toward the center by a factor of 3, apart from a small central 
spike (and corresponding temperature decrease). The initial conditions for the 2D run were 
therefore chosen to be a uniform density with a higher temperature in the jet, in pressure 
equilibrium with the transverse magnetic field changing sign in the center of the jet. The 
code is a 2.5D compressible MHD code which is spectral in one direction (the flow direction 
y) and uses compact finite differences on a non-uniform grid achieving a spatial resolution 
of 0.07 AU in the x direction. The overall number of grid points was 130 in the y direction 
and 128 in the x direction. Non-reflective boundary conditions are imposed in this direction 
using the method of projected characteristics as in Einaudi et al. (2001). Periodic boundary 
conditions are imposed on the y direction. (The x and y axes in the 2.5D run correspond, 
respectively, to the z and x axes in the 3D run.) Initially, the jet-magnetic field equilibrium 
structure is perturbed by random noise. After an initial transient, the jet develops a sinuous 
instability, as may be seen by inspecting the density contours shown for different times in 
Fig. 8. Initially the density is uniform along the jet (the left upper panel). We can see that 
as the time evolves the jet starts to oscillate in a sinuous mode. Because the jet in this case 
is of infinite extent, one can not observe the flag-like oscillation involving a total departure 
from the z = axis as seen in the 3D run. Figure 9 shows the growth rate measured by 
Bettarini (2003). The value obtained is F = 5.0 x 10^ sec""^ for the maximum wavelength 
A = 25 AU . These values are almost identical to what we observe in the 3D runs, for the 
growth rate and for the A, which reinforces the idea that the jet oscillation is due to the 
development of the KH instability. 



5. Conclusions and Discussion 

In the present study we increase and extend the spatial resolution obtained in Opher et 
al. (2003) to study in more details the behavior of the jet. The jet-sheet structure is again 
present and extends for 150 AU beyond the Termination Shock. Turbulent vortices with 
sizes ~ 20 AU are present due to the shear between the jet flow and the surrounding flow 
being pushed aside. The jet is unstable as in Opher et al. (2003). We show that the a de 
Laval nozzle occurs at the TS near the current sheet and it is responsible for accelerating 
the flow to supersonic values after the shock. We address how the width of the current sheet 
depends on spatial resolution showing that as the resolution increases, the width tend to 
a finite value of 4 AU . As we increase the spatial resolution, the profiles of the velocity 
tend to a common shape. This result we believe, indicate that the width of the jet in the 
high resolution calculation, is independent of the grid resolution, and depends on physical 
conditions rather than on numerical resolution. 
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In order to verify and confirm that the instabihty of the current sheet is a KH type of 
instabihty we have also run a 2D compressible MHD case. We used a simple geometry of an 
infinite jet-magnetic field equilibrium using the profiles for the magnetic field and velocity 
from our full 3D runs. The jet is unstable and develop a sinuous instability as observed in the 
3D run. The growth rate is 5.0 x lO^^sec"^ in excellent agreement with the value measured 
in the 3D runs (5.4 x 10~^ sec^^). The growth rate had a maximum for A = 25 AU again in 
agreement with the A measured in the 3D runs, 23 AU. This reinforce the idea that the jet 
oscillation is due to a KH instability. 

The presence of the extended jet, turbulence and the magnetohydrodynamic instabilities 
reinforce the view that the boundary between the solar wind and the interstellar wind is a 
dynamic environment. It is highly distorted discontinuity and that magnetic effects are a 
major player in this region. Still other important effects such as neutrals and the tilt of 
the magnetic axis in respect to the rotation axis need to be included before we can address 
definitively how this boundary appears. Also, when the interstellar magnetic field is included, 
reconnection between the solar and interstellar magnetic field will occur. 

The inclusion of a tilted heliospheric current sheet very likely will introduce qualitative 
changes in the picture represented in this paper. Nerney, Suess & Schmahl (1995) investi- 
gated analytically the solar cycle imprint on the Heliosheath. They predicted that the 11-year 
solar cycle would imprint an alternating magnetic polarity envelopes in the Heliosheath. On 
a much finer scale, there will be strongly mixed polarity regions between the magnetic en- 
velopes due to the reverse magnetic field polarity due to the 25.5-day solar rotation period. 
The velocity shear instability requires hundreds of years to develop, while the solar magnetic 
field changes polarity in much shorter time scales. Thus, in order to have a complete picture 
of the heliosheath, time-dependent boundary conditions on the solar magnetic field need to 
be included. By the picture described in this paper, we expect however, that this region will 
be very turbulent and dynamic. 

The Voyager spacecraft is now at ~ 85 AU and is approaching the Termination Shock 
(there is current controversy on whether it already encounter it (Krimigis et al. 2003; Mc- 
Donald et al. 2003)). Voyager travels 35° above the ecliptic plane towards the nose of the 
heliosphere. We expect that Voyager will encounter an extremely turbulent dynamic region 
that will have a major effect in the energetic particle detected by the instruments on Voyager. 

Finally, for other magnetized rotating stars in motion with respect to the interstellar 
medium we expect to find similar phenomena. We expect the presence of an extended jet 
or disk with length of ~ 150 AU and a width of 10 AU beyond the location where stellar 
winds are shocked in respect to the interstellar medium that surround them. 
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APPENDIX: DEPENDENCE OF THE WIDTH OF THE JET ON SPATIAL 

RESOLUTION 

In this appendix we present the grids for the three cases. Figure 10-12 show the contours 
of the azimuthal magnetic field By for the coarse, refined and super-refined-extended cases. 
The red lines indicate the change in the grid resolution. The numbers indicate the resolution 
in each region, in AU. For the coarse grid, the resolution between the TS and the HP, at the 
HCS, is 6.0 AU. For the refined grid, the resolution between the TS and HP, at the HCS, 
is 1.5 AU. In the super- refined-extended case the resolution between the TS and HP, at the 
HCS, is 0.75 AU. 

Figure 13 shows the magnetic field (solid line) and the velocity profiles (dashed line) for 
the super-refined extended case. Notice that the magnetic field profile is inverted. We can 
see how the flow profile adjust itself to the magnetic one. For the three cases the flow profile 
is broader than the magnetic one. As we increase the resolution, the two profiles almost 
overlap, and the width of the jet is shaped by the HCS width. 
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Fig. 1. — (a) Contours of pressure at t = 1.0 x 10^ sec in the meridional plane (x-z) without 
the solar magnetic field. Same as (a) but with the solar magnetic field. The case presented 
here is the coarse case, with resolution of the grid at the current sheet of ~ 3AU. 
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Fig. 3. — Line plot of the equatorial cut at t = 1.5 x lO^sec of the ratio velocity per sounds 
speed a vs. x. Three cases are shown: super-refined-extended (red), refined (green) and 
coarse (blue). The green and blue refers to the cases at Opher et al. 2003. As we increase 
the resolution the jet(red) extends farther away. The de Laval nozzle is still present. The 
refined jet (green) lose its power at x = —230 AU (indicated by the green arrow), while the 
change in resolution from 1.5 to 3.0 AU occurs at x = —320 AU . The extended jet (red) 
lose its power at x = —320 AU (indicated by the red arrow) while the change in resolution 
(from 0.75 to 1.5 AU ok. x = -350 AU and from 1.5 to 3.0 AU ok. x = -365 AU . 
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Fig. 4. — (a) Schematic drawing illustrating the formation of the de Laval Nozzle, (b) Close 
view of the region of the de Laval nozzle just after the TS. The red arrow indicate the de 
Laval nozzle, where the flow converges. 
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Fig. 5. — (a) Width of the current sheet and the jet vs resolution measured at x 
The red hne is the magnetic field width and the blue line is the velocity width 
(b) Velocity profiles for the coarse, refined and super-refined-extended cases. 
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Fig. 6. — (a) Displacement of the current sheet (AU) vs t/lO^sec for x = —225AU (red) 
and X = —290AU (green). Indexes (I)-(III) marks the different regions of time scales of 
oscillation, (b) Contours of the magnetic field at t = 16.8 years. The scale is the same as 
Fig2a. 
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Fig. 7. — (a) Initial profile for the velocity of the jet VOy as a function of latitude x for the 
2D simulations, (b) Initial profile for the magnetic field BOz as a function of the latitude x. 






Fig. 8. — Density contours for the 2D compressible MHD run in the xy coordinate system. 
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Fig. 9. — Displacement of the current sheet vs t for the 2.5D compressible MHD run. 
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Fig. 10. — The contours of the azimuthal magnetic field, By, for t = 10^ sec for the coarse 
case. The red fines indicate tfie cfianges in grid resolution. Tfie numbers indicate tfie grid 
resolution (in AU). 
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Fig. 11. — The contours of the azimuthal magnetic field, By, for t = 10^ sec for the refined 
case. The red fines indicate the changes in grid resolution. The numbers indicate the grid 
resolution (in AU). 
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Fig. 12. — The contours of the azimuthal magnetic field, By, for t = 10^ sec for the super- 
refined-extended case. The red fines indicate the changes in grid resolution. The numbers 
indicate the grid resolution (in AU). 
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Fig. 13. — Velocity and magnetic field profiles for the super-refined-extended case. The full 
line is the magnetic field and the dashed line is the velocity profile. Note that the magnetic 
field profile is inverted. 



